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Abstract 

We propose a new algorithm for computing validated bounds for the so¬ 
lutions to the first order variational equations associated to ODEs. These 
validated solutions are the kernel of numerics computer-assisted proofs 
in dynamical systems literature. The method uses a high-order Taylor 
method as a predictor step and an implicit method based on the Hermite- 
Obreshkov interpolation as a corrector step. The proposed algorithm is 
an improvement of the C^-Lohner algorithm proposed by Zgliczyriski and 
it provides sharper bounds. 

As an application of the algorithm, we give a computer-assisted proof 
of the existence of an attractor set in the Rossler system, and we show 
that the attractor contains an invariant and uniformly hyperbolic subset 
on which the dynamics is chaotic, that is, conjugated to subshift of finite 
type with positive topological entropy. 
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1 Introduction. 


The aim of this paper is to provide an algorithm that computes validated en¬ 
closures for the solutions to the following set of initial value problems 


r x(t) = /(x(t)), 

I Vit) = Dfixit))-Vit), 

I x(o) e [xo], 

[ ViO) G [Fo], 


( 1 ) 


where / : R” —>■ R” is a smooth function (usually analytic in the domain) and 
[xo] C R", [Vb] C M" are sets of initial conditions. In contrast to standard 
numerical methods, one step of a validated algorithm for © produces sets [xi] 
and [Vi] that guarantee to contain x{h) G [xi] and V{h) G [Vi] for all initial 
conditions x(0) G [xq] and 1^(0) G [Vb], where h > 0 is a time step (usually 
variable) of the method. The computations are performed in interval arithmetics 
in order to obtain guaranteed bounds on the expressions we evaluate. 

The equation for V(t) is called the variational equation associated with an 
ODE. Solutions V{t) give us an information about sensitivities of trajectories 
with respect to initial conditions. They proved to be useful in finding periodic 
solutions, proving their existence and analysis of their stability [3111HKH 
USKE]. They are used to estimate invariant manifolds of periodic orbits [IIS]- 
Derivatives with respect to initial conditions are used to prove the existence 
of connecting orbits IIISIIIII or even (non)uniformly hyperbolic and chaotic 
attractors |40l [41]. First and higher-order derivatives with respect to initial 
conditions can be used to study some bifurcation problems inilillil. This 
wide spectrum of applications is our main motivation for developing an efficient 
algorithm that produces sharp bounds on the solutions to (P). 

In principle, the problem o can be solved by any algorithm capable to 
compute validated solution to IVP for ODEs. There are several available algo¬ 
rithms and their implementations — just to mention a few of them: VNODE-LP 
[371 m [331130] , COSY Infinity P [331133], CAPD P, Valencia-IVP [33]. The 
above mentioned ODE solvers are internally higher-order methods with respect 
to the initial state, which means that they use at least partial information about 
the derivatives with respect to initial conditions to reduce the wrapping effect. 
Therefore, direct application of these solvers uses a higher effective dimension 
(the internal dimension of the solver) than the dimension of the phase space. 
In the case of the codes VNODE-LP and CAPD, this effective dimension is 
{n{n -I- I))^, which dramatically decreases the performance of these methods 
when applied directly to the extended system ©• This key observation mo¬ 
tivated developing the C^-Lohner algorithm [47], which takes into account the 
block structure of © and works in effective dimension. Even in low dimen¬ 
sions, it is orders of magnitude faster than direct application of a solver to 
the variational equations. However, it does not use derivatives of V(t) with 
respect to other coefficients of V (i.e. second order derivatives of the original 
system) to better control the wrapping effect. This is why it usually produces 
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worse estimations than those obtained from direct application of a solver to 
the extended system. 

In this paper, we propose a new algorithm for computation of validated 
solutions to ©• Our algorithm consists of two steps. First, the high-order 
Taylor method is used as a predictor step. Then, an implicit method based 
on the Hermite-Obreshkov (HO) formula is used to compute tighter bounds for 
the variational equations. This last step is motivated by the very famous and 
efficient algorithm proposed by Nedialkov and Jackson |29j and implemented by 
Nedialkov in the VNODE-LP package [55]. We name the proposed algorithm 
C^-HO because it computes bounds for the first order variational equations and 
it is based on the Hermite-Obreshkov interpolation formula. 

Our algorithm, by its construction, cannot produce worse estimations than 
the C^-Lohner algorithm. Complexity analysis (see SectionjSj) shows that, in low 
dimensions, it is slower than the C^-Lohner algorithm by the factor 9/8 only. 
This lack of performance is compensated by a significantly smaller truncation 
error of the method. This allows to take larger time steps when computing 
the trajectories and thus our algorithm appears to be slightly faster than the 
C^-Lohner in real applications — see Section |5| for the case study. 

We would like to emphasize that the proposed method can be directly ex¬ 
tended to the nonautonomous case without increasing the effective dimension 
of the problem. For simplicity in the notation, we consider the autonomous 
case only. In |5], we provide an implementation of the C^-HO algorithm for the 
nonautonomous case. 

As an application of the proposed algorithm we give a computer-assisted 
proof of the following new result concerning the Rdssler system [36]. 

Theorem 1 For the parameter values a = 5.7 and b = 0.2 the system 

{ X = -y- z 

y = x + by (2) 

i = b + z{x — a) 

admits a compact, connected invariant set A that is an attractor. There is an 
invariant subset TL G A on which the dynamics is uniformly hyperbolic and 
chaotic, that is, conjugated to a subshift of finite type with positive topological 
entropy. 

Verification that an ODE is chaotic is not an easy task in general. After de¬ 
velopment of rigorous ODE solvers there appeared numerous computer-assisted 
proofs of the existence of chaos in classical low-dimensional systems — just to 
mention two pioneering results [24l [46]. To the best of our knowledge there 
are only two computer-assisted proofs of the existence of chaotic and (non)- 
uniformly hyperbolic attractors for ODEs gniEi]- These results became pos¬ 
sible with development of suitable theory and the algorithms capable to inte¬ 
grate variational equations. In |41j the C^-Lohner algorithm implemented in the 
CAPD library [6] was used. 

The paper is organized as follows. In Section ITtI we introduce the symbols 
and notation used in the paper. Section [5] contains description of the algorithm 
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and the proof of its correctness. In Section [31 we analyze the complexity of 
the C^-HO algorithm and we compare it to the complexity of the C^-Lohner 
algorithm. In Section m we compare the bounds obtained by the C^-Lohner and 
C^-HO algorithms on several examples. In Section [SI we give a more detailed 
statement and proof of Theorem |T] We discuss also how the computing time 
depends on the choice of the C^-Lohner and C^-HO algorithms to integrate 
variational equations. 

1.1 Notation. 

By I we denote the identity matrix of the dimension clear from the context. By 
Df we denote the derivative of a smooth function / : R" —>■ R™. By D^f we 
denote the partial derivative of / with respect to x. 

The local flow induced by an ordinary differential equation (ODE) x{t) = 
f{x{t)) will be denoted by ip, i.e. (p{-,x) = x(-) is the unique solution passing 
through X at time zero. We will often identify an element x G R" with the 
function x(-) = (p(-, x). We denote x, V) = Dx<p{t, x) ■ V. Clearly '0(-, x, V) 
is a solution to the first-order variational equation associated with an ODE with 
the initial conditions x and V. 

Let / : R ^ R” be a smooth function. By ( x) we denote the vector of fth 
derivatives of /. Normalized derivatives (Taylor coefficients) will be denoted by 
/W(a:) = i/W (a;). We apply this notation to derivatives and Taylor coefficients 
of the flows if and tp taken with respect to the time variable 

:= ((/?W (•, a;))(t), 

E) := x,V)){t). 

Interval objects will be always denoted in square brackets, for instance [a] = 
[a, a] is an interval, and [u] = ([ui],..., [un]) is an interval vector. Matrices or 
interval matrices will be denoted by capital letters, for example [A]. Vectors 
and scalars will be always denoted by small letters. We also identify Cartesian 
product of intervals [ui] x • • • x [u„] with a vector of intervals ([ui],..., [u„]). 
Thus interval vectors can be seen as subsets of R”. The same identification will 
apply to interval matrices. 

The midpoint of an interval [a] = [a, a] will be denoted by a = (a-|-a)/2. The 
same convention will be used to denote the midpoint of an interval vector or 
an interval matrix; for example for [u] = ([ui],..., [?;„]) we put v = {vi,... ,Vn)- 
Sometimes, we will denote the midpoint of products of interval objects by 
mid([E][r]). 

Throughout this article, interval vectors or interval matrices marked with 
tilde [^, [V] will always refer to rough enclosures for the solutions to the IVP 
problem o — see Section [2T] for details. 
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2 The algorithm. 

Consider the initial value problem o and assume that we have already proved 
the existence of the solutions at time tfe, and we have computed rigorous bounds 
[xk] and [14] for ip{tk,[xo]) C [xk], 'ip(tk,[xo],[Vo]) C [Vk], respectively. Let us 
fix a time step hk > 0. A rigorous numerical method for o consists usually of 
the following two steps: 

• computation of a rough enclosure. In this step, the algorithm validates 
that solutions indeed exist over the time interval [4,4 + and it pro¬ 
duces sets [^ and [V], called rough enclosures, which satisfy 

(p{[0,hk],[xk]) C [y] and (3) 

i;{[Q,hk],[xk],I) C [Cj. (4) 

• computation of tighter bounds [xfc+i], [14-i-i] satisfying (^(4 + hk, [a;o]) C 
[xk+i] and i>{tk + hk, [xo], [lb]) C [Vfc+i]- 

In the sequel, we give details of each part of the proposed algorithm. 

2.1 Computation of a rough enclosure. 

This section is devoted to describe a method for finding rough enclosures (H 
0]). The key observation is that the equation for V in ([1]) is linear in V, which 
implies that the following identity holds 

^/’(4 + hk,Xo,Vo) = tfj{hk,(p{tk,xo),I) ■ 'ilj{tk,Xo,Vo) 

provided all quantities are well defined. This implies that 

f/’(4 + hk, [xo], [Vb]) C tpiyhk, [xk\,I) ■ [Vk\. (5) 

Hence, it is sufficient to use I as an injtial condition for the variational equations 
when computing a rough enclosure [V]. 

One approach to find rough enclosures [^ and [H] is to compute them sep¬ 
arately. Given a set [^ satisfying (jS]) and computed by any algorithm [26l [30] , 
we can try to find an interval matrix [V] such that 

I+[0,hk]Dfm-[V]c[V]. (6) 

If we succeed, then the interval matrix [H] satisfies (|4|). This method is known 
as the First Order Enclosure (FOE). It has, however, at least one significant 
disadvantage. If we already know that the solutions to the main equations exist 
over the time step hk ([^ has been computed) there is no reason to shorten this 
time step because solutions to variational equation also exist over the same time 
range. However, this shortening might be necessary to fulfill the condition (|5)) . 
To avoid this drawback, Zgliczyhski proposes |47] a method based on logarithmic 
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norms that always computes an enclosure [y] satisfying (|3]) for the same time 
step hk, provided we are able to find an enclosure satisfying ([3]). This type of 
enclosure is also used in the C’’-Lohner algorithm |45] for higher order variational 
equations. 

Another strategy is to use the High Order Enclosure (HOE) method [TOl HZl 
ISO]. The authors propose to predict a rough enclosure of the form 

m 

[^=^[0,MV^’(0,[^.]) + H, (7) 

2=0 


where [s] is an interval vector centered at zero. The inclusion 


[0,hfc]-+V'”^+''(0, [^)c[e] 


( 8 ) 


implies that the set is indeed a rough enclosure, i.e. it satisfies ([3]). If the 
inclusion ([5]) is not satisfied, then we can always find hk < hk such that ([5]) 
holds with this time step, and thus is a rough enclosure for the time step 
hk- This strategy is very efficient because we do not need to recompute 
Furthermore, with quite high order m and a reasonable algorithm for time step 

l/(m+l) 

« 1 . 


prediction, we usually have hk/hk= ||y|m+i!(o,[y])|| ) 

The above method can be used to find simultaneously two enclosures I 
for the entire system ©• We predict [j/] as in ([7]) and 


],[r]) 


m 

[r] = E [O’/) + [£;]. (9) 

2=0 

Then we check simultaneously (HJ and 

[0,M-+V[-+il(0,[^,/)[H] C [E]. (10) 


Due to linearity of the equation for variational equations we can consider two 
strategies when OT is not satisfied. We can 

1. either shorten the time step as we do in (|5]) or 

2. compute [H®] from the logarithmic norms for the same time step hk as 

proposed in [47] and set [E] = [0, (0, [^,/)[H°]. Then [V] 

computed as in (jH]) with such [E] satisfies (|3|). 

We would like to emphasize that in both cases we do not need to recompute the 
coefficients '0^(0, [^,1) which is very expensive. The first strategy is recom¬ 
mended when the tolerance per one step is specified which means that there is 
a maximal norm of [E] which should not be exceeded. The second strategy ap¬ 
plies when the fixed time step is used (by the user choice or application specific 
reason). 
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Algorithm 1: Predictor. 

Input : m - natural number (order of the Taylor method) 
hk - positive real number (current time step) 

[xk], - interval vectors 

[V] - interval matrix 

Output: ([ccg+J, [r°], [P°], [R°]) C K” x R" x R"" x R”" 

Compute: 

yo ^ T,Zo 

[y]^Etohl^^^iod^k]y, 

[r°] ^ [0,/ife]'"+Vt'"+^'(0,[^); 

K+i] ^ (yo + - Xk)) n[y] + [r°]; 

[PO] ^ [A] + [R% 

return ([x^+J, [r°], [P°], [i?°]); 


Our tests show that the version of (HOE) gives better results than the 
approach proposed by Zgliczyhski, which uses logarithmic norms. Since com¬ 
puting a rough enclosure 0 is not the main topic of the paper we omit details 
here. 

In what follows we assume that we have a routine that returns three quan¬ 
tities: hk, [H] for which the properties ([3]) and (j4]) hold. 

2.2 The predictor step. 

We give a short description of the C^-Lohner algorithm m which will be used 
as a predictor step in the C^-HO algorithm. 

Lemma 2 Assume hk, [xk], [y\, [V] axe such that |3) and ^ hold. Then the 
quantities ([a:°_|_j^], [r°], [H°], [t?°]) computed by the Algorithm\li satisfy 


(p{hk,[xk]) C [xg+i], (11) 

ib{hk,[xk],I) C [H°], (12) 

[0,M™+V‘"^+'’(0,[^) C [rO] and (13) 

C [R% (14) 

Proof: The Taylor theorem with Lagrange remainder implies that for all Xk S 
[xk] and each component j = 1,..., n: 

m 

(Pj{hk,Xk) = '^hlipj^^yO,Xk) + hf+^ipj^"^+^yTj,Xk) (15) 
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for some Tj G {0,hk). By the assumptions (p(Tj,Xk) G and by the group 
property of the flow, we have 




Therefore 

m 

^{h,,x,) e + hrV'™+''(0,[^) c [y] + [r°]. 

i =0 

Since [xk\ is convex, we can apply the mean value form to the polynomial 
part of m and obtain that for Xk G [x^] there holds 

m m 

hl(p^^^(0,Xk) G hl.(p^^^(0,Xk) + [A]{xk - Xk) C yo + [^]([a;fc] - ®fc)- 

i —0 2—0 

Gathering the above together, we obtain 

^{hk, [Xfc]) C (yo + [^]([a;fc] - %)) G [y] + [r°] = [x^+J. (16) 

In a similar way, we deduce that for Xk G [x^] and for each component j, c = 
1,..., n there holds 

G [^, [F]) 

for j, c = 1,..., n, and in consequence 

i;{hk,[^k],I)c[A] + [R^] = [Vy 


2.3 The corrector step. 

The goal of this section is to set forth a one-step method that refines the results 
obtained from the predictor step and returns tighter rigorous bounds for the 
solution to the ODE and its associated variational equation (ED- The method 
combines the algorithm by Nedialkov and Jackson [5U] based on the Hermite- 
Obreshkov interpolation formula with the C^-Lohner algorithm for variational 
equations proposed by Zgliczyhski m- For reader’s convenience, we recall here 
the key ideas of the Hermite-Obreshkov method. 

For natural numbers p, q, i such that i < q, let 


„9.P 



For a smooth function u : K R" and real numbers h,t we define 

'^q^p{h,u,t) = (t). 

i=0 
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Using this notation, the Hermite-Obreshkov [32] formula reads 

4',,p(-h, u, h) = u, 0) + {-lYc^/hP+<i+^R{h, u), (17) 

where 

R{h,u) = , n € {0,h),i = 

The key observation which was the main motivation to develop rigorous 
numerical method based on this formula is that the coefficient 
can be very small for p = q. Thus, this formula can have significantly smaller 
remainder than the Lagrange remainder used in the Taylor series method. 

Now we would like to apply (inD to the flows ip and ip := Let [xk\ 

be a set of initial conditions and assume that from the predictor step we have 
computed ([a;^+i], [r°], [U°], satisfying (fTTMTil) . 

Let us fix positive integers p, q such that m = p + q, Xk G [a;^] and put 
Xk+i = ip{h,Xk)- The formula (flTl) applied to this case reads 

J24’’'i-hky<P^"H0,Xk+i) = +e, 

i—0 z—0 

where e G (—Identifying vectors Xk, Xk+i with unique solutions 
Xk{-), Xk+i{') to the ODE passing through them at time zero, we obtain the 
equivalent but shorter form 

^q,p{-hk,Xk+i,0) = ^p^q{hk,Xk,0) + e. (18) 

Take the midpoints G [a^fc+iJi G [a;^]. Since interval vectors are 

convex sets, and the local flow is a smooth function in both variables, we can 
apply the mean-value form to both sides of (fTSl) to obtain 

lllg,p( Xk-^-iT “t“ J— (a^fc-t-1 a^/g_|_i) — p^q{hk ; a^fc; 0) T J-\- (a^fc a^fc) T ^ 

for some 

J_ G [D^^q^p{-hkyxl^y\,Q)\ and 

J+ G \Dx^p,q{]Xk^\Xk\i^y\ • 

We obtained a linear equation for Xk+i 

a^fc_|_i) — J-\-{Xk a^fc) T (Ikp^g(/i/j;, Xfc, 0) ^q,p(, 0))T^ (^^) 

in which the matrices J± are unknown, but they can be rigorously bounded. 
Denoting 

[( 5 ] = '^p^q{hk,Xk,0) — 'i^q,p{—hk,Xk^i,0), 

[e] = (-l)fo?’"[r°], 

[•/-] = [Dx'i'q,p{—hk,[x^_^i],0)\, 

['^+] ~ P^qi^kj [xk]j^)] ^ 

[5] = i-jzyj-], 

H = JiMH + N) + [5]([x°+i]-s°+i) 
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and applying the interval Krawczyk operator [IJII1I3I] to the linear system 
mi), we obtain that for Xk S [a;^], 

ip{hk,Xk) = Xk+i G (Nfe] - ®fc) + [r] (20) 

which is the main evaluation formula in the interval Hermite-Obreshkov method 
for IVPs presented in [29] . Note that this formula has exactly the same structure 
as used in the predictor step. 

Each coefficients of [S'] is an interval containing zero and its diameter tends 
to zero with hk —>■ 0. The vector [^] is almost a point vector, and [e] can be 
made as small as we need (manipulating the time step hk). Therefore, the total 
error accumulated in [r] is usually very thin in comparison to the size of the set 
[xk] we propagate. Thus, the main source of overestimation when evaluating 

(HOD comes from the propagation of the product ^J_^[J+]^ ([cc/c] — Xk). There 
is a wide literature on how to reduce this wrapping effect for such propagation 
(see [26] for a survey), and we will give some details concerning this issue in 
Section 12.41 

In what follows we argue that, with a little additional cost, we can compute 
a possibly tighter enclosure for the solutions to variational equation than the 
bound [E°] obtained from the predictor step. Let us fix Xk G [xk] and let 

V = ipihk, Xk, I). Applying (I17|) to the solutions to the variational equation, we 
obtain 

^ cf P(-hfc)V'*’ (0, Xk+i,V)=Y, (0, Xk,I) + E, 

i—0 z—0 

where E G (—l)‘^c^’^[i?°]. Since ip is linear in V, we obtain that the matrix 

V = ip{hk, Xk,I) belongs to the solution set to the linear equation 

[J.]V = [J+] + [E], 

where [E] = (—l)'?c^’P[i?°]. Note that from the predictor step we already know 
that V G [E°]. Applying the interval Krawczyk operator [1] [181 [31] to this linear 
system we obtain 

K G + [E]) + (1- J-V-mV°] = + [E]) + [5][K°]. 

Due to linearity of the variational equation, we can reuse the matrices [J-], [^+], 
[S'] and JZ^ computed in the corrector step for ip. Thus the additional cost is 
just a few matrix additions and multiplications. Algorithm [5] and Lemma [3] 
summarize the above considerations. 

Lemma 3 Assume that hk, [xk], [^, [K] are such that and ^ hold and 
that the quadruple ([a;^_|_]^], [r°], [K°], [A°]) is returned by the predictor step (Al- 
gorithm\^. Then the quantities ([xfe+i],[K]) computed by Algorithm\^ satisfy 

p{hk,[xk\) C [Xk+i] and (21) 

ip{hk,[xk],I) C [K]. (22) 
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Algorithm 2: Corrector. 

Input ■ p,q - positive integers 

hk - positive real number 
[xk] - interval vectors 

Q^fc+i]> [^°]> [^°]) [^°]) " from the predictor step with 
m = p -\- q 

Output: ([xfe+i], [F]) 

Compute: 

fel 

[JJ\^[Dj,,p{-hkAxl+^],Q)\-, 

[d+] ■<— [Dx'^p,q(hk,{xk\,0)]] 

[S]^I-JZ\J-\, 

[r]^JzH[5] + [e]) + [S]{[xl^,]-^k+i)-^ 

[i?]^ Jri([J+] + (-l)«c«'P[i?0]); 

[xk+i] ^ (x^+i + (jr^[J+]) ([xfc] - Xk) + [r]) n 
[F] ^ {[R] + [5][C°]) n [C°]; 

return ([x^+i], [C]); 


We would like to emphasize that by its construction the proposed algorithm 
always returns tighter bounds than the C^-Lohner algorithm because the result 
obtained from the corrector step is intersected with the bound obtained from 
the predictor step. 

2.4 Propagation of product of interval objects. 

It is well known that evaluation of expressions in interval arithmetic can produce 
large overestimation due to dependency of variables and the wrapping effect 
[2 HOI |25l [3l]. To reduce this undesirable drawback we follow the ideas from 
[20l [26l [29l |47] , and we represent subsets of M" and M" in the forms (doubletons 
in [26] terminology) 


[xfc] — Xk -t- [sfc] and (^1^) 

[Vk] = Vk + Ak[Rk] + Qk[Sk]. (24) 

The initial conditions ([xq], [Vq]) of ([T]) are assumed to be already in the form 
((23H24)) . The parallelepipeds Xk + Ck[rk] and 14 + Ak[Rk] are used to store 
the main part of the sets [x^] and [Vfc], respectively. The terms Bk{sk\ and 
Qk[Sk\ are used to collect all usually thin quantities that appear during the 
computation. 

According to ([5]), the bound for ip{tk + hk, [xq], [Vq]) can be computed as 

[Pfe+i] = [C][14], 
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where [F] satisfies (1^ . Substituting the representation (1^ we obtain 
[^fe+i] c [y] (14 + ^fc[-Rfc] + Qk[Sk]) n (Vfc+i + Ak+i[Rk+i] + Qfc+i[iS'fc+i]), 
where the new representation is computed as follows 
[AA] = (^[V]-v){Vk + Ak[Rk]), 

14+1 = VVk, 

Ak+i = V Ak, 

[^fc+i] = and 

[^fc+l] = [Rk]- 

In principle, the matrix Qk+i can be chosen as any invertible matrix. The nu¬ 
merical experiments [iniissiiiniiiT] show that one of the most efficient strategies 
in reducing the wrapping effect is to compute Qk+i as an orthogonal matrix 
from the QR decomposition of the point matrix VQk- Note, that even if the 
matrix Qk+i is a point matrix, the inverse Qk^i must be computed rigorously 
in interval arithmetic. 

Similar strategy is used for propagation of products in 
[Xk+i] c -f (jr^[J+]) {[Xk] - Xk) + [r] 

= {Ck[rk\ + Bk[sk]) + H 

— see [201 EH EH Hz] for details. 

3 Complexity. 

In this section, we explain why the C^-HO algorithm may perform better than 
the C^-Lohner algorithm, even if it has higher computational complexity. A 
large numerical and theoretical study were performed to compare the Interval 
Hermite-Obreshkov method (IHO) with the Interval Taylor Series Method (ITS) 
m- It has been shown that, with the same step size and order, the IHO method 
is more stable and produces smaller enclosures than the ITS method on constant 
coefficient problems. Furthermore, the IHO method allows the use of a much 
larger stepsize than the ITS method, thus saving computation time during the 
whole integration. However, comparing these two methods in the nonlinear 
case is not as simple as in the constant coefficient case. Our goal is to predict 
the benefits of performing additional calculations required by the IHO method 
applied to O- 

3.1 Cost of C^-Lohner and C^-HO methods per step. 

We assume that both predictor (Algorithm [1]) and corrector (Algorithm [2]) have 
the same order. That is, if the order of the predictor is m, we consider the 
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corrector step with p and q such that m = p + q. In what follows we list the 
most time-consuming items of the predictor and corrector, which are the core 
of their computational complexity. 

In the analysis give below, we count the number of operations which are 
really executed by the implementation, rather than the possible theoretical and 
asymptotic complexity. Therefore, we assume that the product of two square 
interval matrices is computed by the naive algorithm (three nested loops or 
equivalent), which executes exactly interval multiplications. 

We would like to emphasize, that the rigorous integration of a differential 
equation is a very difficult task even in quite low dimensions. Thus, dimensions 
used in practice are usually less than 20. Computer-assisted proofs for 100- 
dimensional systems are actually the state of the art — see for instance [15] . 
Therefore, the use of asymptotically fast algorithms for matrix multiplications, 
such as the Strassen algorithm [35] or the Coppersmith-Winograd [3] algorithm, 
does not make any sense. 

Let us denote by c/ the cost of evaluating the vector field ([T|). For the C^- 
Lohner step we need the following operations (predictor step and propagation 
of doubleton representations) 

• simultaneous computation of (0, [xs,]) and (0, [xk],I) up to order m. 
This is performed by means of automatic differentiation techniques, and 
it takes c/(2n -|- l)(m -I- l)(m -I- 2)/2 multiplications — see [53] . 

• simultaneous computation of (0, [y\) and (0, [y],I) up to order m-l-l. 
This is performed by means of the automatic differentiation techniques and 
it takes c/(2n -|- l)(m -|- 2)(m -|- 3)/2 multiplications — see [53] . 

• 13 matrix by matrix multiplications, 2 point matrix inversions and 2 point 
matrix QR decompositions. Approximate QR decomposition of a point 
matrix is much cheaper than the product of interval matrices and we may 
assume that it takes O(n^) with a constant less than one (in terms of in¬ 
terval multiplications). The inversion of a point matrix which is very close 
to orthogonal is performed by means of the interval Krawczyk operator 
[l8] and takes (one multiplication) because an approximate result is al¬ 
ready known (transposition of an approximate orthogonal matrix). Thus, 
the total cost of all matrix operations listed above is at most 17n^. 

We did not list cheaper operations like additions, intersections of interval ob¬ 
jects, matrix by vector products. All polynomial evaluations perform in total 
0{n?m) interval multiplications, and they add significant cost to linear systems 
(c/ = 0) and to nonlinear systems but with very small number of nonlinear 
terms (c/ -C n). Thus, we skip them. 

To sum up, the total costs of the C^-Lohner step is 

C-Loin, m) ~ Cf{2n -I- 1)(to -I- 2)^ -|- 17n^. 

In the C^-HO method, we can reuse the Taylor coefficients of (f and ip com¬ 
puted in the predictor step, which are needed for computing the [J+] matrix 
and 'I'p,q(/ifc, Xfc, 0). Thus, the additional cost is 
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• computation of [s^fc+i]) I) up to order q. This is performed by means 

of automatic differentiation techniques, and it takes c/(2n + 1 )(( 3 ' + 1 )(( 3 ' + 
2)/2 operations — see [33], 

• computation of /ifc,0) takes Cf{q + 1)(9 + 2)/2, 

• rigorous inversion of the point matrix J_ takes at most 2n^ (one non¬ 
rigorous inverse and one interval matrix multiplication in the Krawczyk 
method) and 

• 4 interval matrix multiplications require in total 4n^ operations. 

The total additional cost of the C^-HO step is at most 

Cf{n + l){q + l){q -f 2) -|- 6n^ . 




Figure 1: Plot of C}io{n,m)/Ci,o{n,rn) for Cf =n and c/ = respectively. 

Assume that m is an even number and take q = p = ^. Then the above 
additional cost of the C^-HO method is approximately 

ic/(n -I- l)(m -I- 2)(m -f 4) -|- 6n^. 

Hence, total computational complexity of the C^-HO method is 

CYio{n,m) cs Ci,o{n,m) + ^c/(n -|- l)(m -f 2)(m -I- 4) -|- 6n^. 

In general, the complexity depends on the cost of the vector field evaluation c/ 
which can be arbitrary. In Fig. [T] we plot the graph of Cho/Clo for two cases. 
The case c/ = n means that the number of nonlinear terms in the vector field 
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is equal to the dimension of the problem. We observe that, if order m of the 
method is much smaller than the dimension n, then the complexity is dominated 
by the matrix operations and we have 
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lim Ciioi'n,m)/Ci^oin,m) = — ~ 1.35294. 

n—^oo 17 

for all fixed values of m. We observe, however, that for reasonable dimensions 
and orders, this factor is much smaller than the limit value. 

A model example for the Cf = n? case is a second order polynomial vector 
field with nonzero coefficients in the quadratic terms. In this case we have 


lim CHo{n,m)/CLo{n,m) 

n—>oo 


9m^ + 38m + 132 
8m^ + 3m + 100 


The above analysis shows that the additional cost of the C^-HO method in 
a typical nonlinear case approaches 1/8. In the next section, we argue that this 
extra cost of the C^-HO method is compensated by the larger time steps this 
method can perform without losing the accuracy. 


3.2 Maximal allowed time step for a fixed error tolerance. 

To obtain insights into the compared methods, we ask the following question: 
given an acceptable tolerance e per step, what is the maximal time step h of both 
methods that guarantees achieving this constraint. For the C^-Lohner method, 
we have to solve the following inequality 


/i™+V[™+i](0, [y]) 


< e. 


In general, it is very difficult to answer this question because = [y{h)] de¬ 
pends on h. If [xk\ is a point and e is very small, we can assume that the vector 
field is almost constant near [xk] and thus [^) « (^[’"+^1(0, [xfc]). Since 

[xk] C [^, we always have [a;fc])|| < [y])||. With this sim¬ 

plification, we obtain an upper bound for the time step 


/iLO := h = »"+ 


^l^+i](0,[xk]) 


For the C^-HO method we obtain the following upper bound for the time 


step 


^HO '■= h = 


where by |"m/2] we denote the smallest integer not smaller than m/2. Denote 


g{m) := Iiho/^lo 



(25) 
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It is easy to show that 

lim g{m) = 2. 

m—>-130 

In Fig. [21 we observe that the values of g{Tn) rapidly grow for small values of m. 
This is important from practical point of view — even for small order m = 6 the 
C^-HO method allows up to 53% larger time steps than the C^-Lohner method. 
For m = 16 this is 74%. For larger values of the tolerance e, the computed 
enclosure for h = hno is usually significantly smaller than that computed for 
h = IiHOj which affects the norm [^)||. Therefore, the value g{m) is 

a theoretical upper bound for the possible growth ratio of the time step in the 
C^-HO method achievable when e —>■ 0. 

g(m) 

2.5 



1.0 • 


0.5 


10 


15 


20 


25 


30 


Figure 2: Plot of the theoretical maximal factor of maximal time step in the 
C^-HO and C^-Lohner methods for a fixed tolerance — see [2^ 


4 Benchmarks. 

In this section, we present the results of a comparison of the C^-Lohner algorithm 
and the C^-HO algorithm. The structure of the tests is as follows. For a given 
ODE 

• we take an initial condition u which is an approximate periodic orbit for 
the system; 

• we integrate the variational equations along this periodic orbit using the 
C^-Lohner and C^-HO algorithms with the same algorithm for rough en¬ 
closure (HOE), the same order m = p + q of the methods and a constant 
time step h; 

• we compare the logarithm of the maximal diameter of the interval matrix 
[14] (diameter of the widest component) computed by means of the two 
algorithms; and 
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• we repeat the above two steps six times: for two different orders of the 
numerical methods each for three different time steps. 

Fixing the time steps allows us to compare the size of the enclosures returned by 
the two algorithms over the same time step. This will allow us to conclude that 
the C^-HO algorithm can take larger time steps than the C^-Lohner algorithm 
without significant lost of accuracy. The comparison of the two algorithms with 
variable time steps will be given in Section [S] 

The above test is performed for four ODEs: the Lorenz [21] system, the 
Henon-Heiles system the Planar Circular Restricted Three Body Problem 
(PCR3BP), and a 10-dimensional moderately stiff ODE. Below we give initial 
conditions and discuss obtained results. 

The Lorenz system [5T] for “classical” parameters is given by 


{ X = lQ{y-x), 
y = x{28-z)-y, 

z = xy- § 2 . 

The Henon-Heiles system m is a hamiltonian ODE given by 
( X = -x{l + 2y), 

\ y = - 2/(1 + 2 /)- 


(26) 


(27) 


The PCR3BP is a mathematical model that describes motion of a small 
body with negligible mass in the gravitational influence of two big bodies. The 
motion is restricted to the plane, and the two main primaries rotate around 
their common mass centre. The equations for motion of the small body is then 
given by 

ix-2y = D^n{x,y), 

]^y + 2x = Dyn{x,y), 

where 


^ix,y) = 


+ y^) + 


1- y 


y/{x + y)'^ + y^ 


_ y _ 

{x-l + + y^' 


The parameter y, stands for the relative mass of the two main bodies. Eor our 
tests we fixed y = 0.0009537, which corresponds to the Sun-Jupiter system. 

The last ODE is the Galerkin projection of the following infinite dimensional 
ODE 

k—1 CXD 

(1 k ^ ^ ‘2k ^ ^ ^n^n-\-k 

n—1 n—1 

onto (oi,... ,aio) variables. The above system describes solutions to the one¬ 
dimensional Kuramoto-Sivashinsky PDE [I1I33 under periodic and odd bound¬ 
ary conditions, see |48l [49] for derivation. 
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Figure 3: Approximate periodic orbits for (a) the Lorenz system (l26)l . (b) 
Henon-Heiles hamiltonian (P71) . (c) the PCR3BP (1281) and (d) a the 10- 
dimensional Galerkin projection of the Kuramoto-Sivashinsky equation (|29l) . 
respectively. 


We have chosen initial conditions that are close to periodic orbits of these 
systems (see Fig. ED 


= (-2.1473681756955529387,2.078047612582596404,27), 


^Henon-Heiles 


(0.0,0.10903,, 0.5677233993382853,0.0), 


MPCR 3 BP = (0.92080349132074,0.0,0.0,0.1044476727069111) and 

0.2012106 

1.2899797585174486 

0.2012106 

-0.37786628185377774 
-0.042309451521292417 
“ 0.043161614695331821 

0.0069402112803455653 
-0.0041564870501656455 
-0.00079448972725675504 
0.00033160609117820303 


For the two Hamiltonian systems, the coordinates are given in the order (x, y, x, y) 
The orbit mpcrsbp is the well known Li Lyapunov orbit for the Sun-Jupiter- 
Oterma system. In [JS], a computer assisted proof of the existence of a periodic 
solution for the full infinite dimensional system (1291) is given. The projection of 
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this periodic orbit onto the first 10 coordinates is very close to the point uks- 
In fact, due to very strong dissipation, the variables with high indexes have very 
small impact on the dynamics of (1291) . The system becomes very stiff even for 
relative small dimension of the Galerkin projection. 


m =10, h = 0.01 

O h. 

-2 

-4, 

- 6 ^ 



0.0 0.5 1.0 1.5 


m =10, h = 0.03 



m =10, h = 0.05 



m =15, h = 0.01 





Ci-LO 

C-HO 


Figure 4: The results of the tests for the Lorenz system. We plot S{t) = 
logj^g diam([Id(t)]) along trajectory of the point Muorenz integrated with order m 
and with fixed time step h. 

In Figs. |4]-[3 we present results of our numerical experiments. On these 
figures we show plot of S{t) = log^g diam([y(t)]) along an approximate periodic 
trajectory, where diam([V(<)]) is the largest width of coefficient in the interval 
matrix [Id(t)]. We can see that in each case the C^-HO method does not return 
worse results than the C^-Lohner algorithm. This is due to its construction, 
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m =10, h = 0.01 


m =15, h = 0.01 



m =10, h = 0.04 



m =10, h = 0.07 




= 15, h = 0.04 



m =15, h = 0.07 



Ci-LO 

C-HO 


Figure 5: The results of the tests for the Henon-Heiles system. We plot 
S{t) = logiQ diam([F(t)]) along trajectory of the point unenon-Heiies integrated 
with order m and with fixed time step h. 


20 












m =10, h = 0.01 


m =15, h = 0.01 
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Figure 6: The results of the tests for the PCR3BP. We plot S(t) = 
logj^g diam([R(t)]) along trajectory of the point mpcrsbp integrated with order 
m and with fixed time step h. 
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m =10, h = 0.0005 m =15, h = 0.0005 








Ci-LO 

C-HO 


Figure 7: The results of the tests for the 10-dimensional Galerkin projection of 
the Kuramoto-Sivashinsky equation. We plot S{t) = log^g diam([l/(t)]) along 
trajectory of the point uks integrated with order m and with fixed time step h. 
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because the bounds computed in the corrector step are intersected with the 
estimates from the predictor step, which is used in the C^-Lohner algorithm. 
Indeed, in the Algorithm [5] we have 

[F] ^ {[R] + [5][F°]) n [F°]. 

Looking at the columns of Figs. IH [5] and [6l we observe that in each case the 
advantage of the C^-HO method increases when the time step is enlarged, and 
the obtained bounds can be orders of magnitude tighter. This is due to the fact 
that the C^-HO method has times tighter truncation error than the Taylor 
method. To give some numbers, let us take m = 20 which is a typical order 
used in computations. Then p = q = 10 and ~ 5.4 ■ 10“®. 

Increasing the order of the method makes the truncation error of the C^- 
Lohner method smaller when the time step is fixed. Therefore, in the right 
columns of each figure we observe that the C^-Lohner method performs much 
better than in the left column. The C^-HO method, however, still returns tighter 
enclosures and is capable to take even larger time steps without significant lost 
of accuracy. This is especially important for stiff problems, where the time 
steps used by a nonstiff solver cannot be large and thus integration over large 
time interval is very expensive. We would like to emphasize, that the maximal 
possible time step that a rigorous ODE solver can take is limited mainly by 
the possibility of finding a rough enclosure over the time step. To the best of 
our knowledge, the HOE algorithm [3^, which is nonstiff, is one of the most 
efficient. Therefore construction of a general rigorous stiff ODE solver without 
extra knowledge of the system is a challenge. Remarkable exceptions are solvers 
for infinite-dimensional strongly dissipative systems |111148) . where the structure 
of the system is used to construct a dedicated so-called dissipative enclosure. 

In Fig. Owe can see that the C^-HO method can perform much larger time 
steps keeping very good accuracy of computed bounds. 

The bounds obtained by the C^-HO method are tighter than those returned 
by the C^-Lohner algorithm, but as we observed in Section^ the C^-HO method 
is computationally more expensive. In the next section, we argue that this extra 
cost per step is compensated by the larger time steps we can take. 

5 Applications. 

In this section, we present an application of the proposed algorithm to a computer- 
assisted proof of a new result concerning the Rossler system [36]. We focus on 
the comparison of the time of computation needed to prove this result, when the 
C^-Lohner algorithm and the C^-HO algorithm are used to integrate variational 
equations, which are necessary to prove this theorem. 

The classical Rossler system [36] is given by ©■ When the two parameters 
a, b vary, the system exhibits wide spectrum of bifurcations. This system admits 
period doubling bifurcations [44], which lead to chaotic dynamics [46]. In [33] . 
the existence of two periodic orbits was proved by means of the Conley index 
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theory. All the above results about the system (I2|) are computer assisted and 
use rigorous ODE solvers. 


20 

10 

z 

0 


Figure 8: Typical chaotic trajectory of the system ([2]) and a slice of the Poincare 
section 11. 

Let n = {(a;, j/, z) G : a: = 0 and i > 0} be a Poincare section (see Fig. [8]) 
and let P : n n be the Poincare map. Since the x coordinate is equal to zero 
on n, we use only two coordinates (y, z) to describe points on 11. 

Theorem 4 Let Ib = —10.7, xb = —2.3, Im = —8.4, tm = —7.6, In = —5.7, 
rjv = —4.6, Z = [0.028,0.034] and let 

B = [Ib, tb] X Z, 

M = [Im^xm] X Z and 
N = [Zat, tn] X Z. 

For the classical parameter values a = 0.2, b = 5.7 the following statements 
hold. 

• The system ^ admits an attractor. The set B is a trapping region for the 

Poincare map, i.e. P is well defined on B and P{B) C B. In particular, 
there exists a maximal invariant set A = Un>o for the map P that 

is compact and connected. 

• The maximal invariant set for P^ in fVUM, denoted by TL = inv(P^,7VU 
M) G A, is uniformly hyperbolic; in particular it is robust under pertur¬ 
bations of the system. The dynamics of P^ on Ti is chaotic in the sense 
that P^lu is conjugated to the Bernoulli shift on two symbols. 

Proof: The tools used in a computer-assisted proof of Theorem |4] are well 
known, and we summarize them here. 
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Figure 9: Geometric conditions that guarantee the existence of chaotic dynamics 
— see (1^ . 


Trapping region. Verification that B is a trapping region for P reduces to 
checking the inclusion 

PiB) C B. 

We uniformly subdivided the set B onto N = 160 pieces of the form Bi = 
[Ui, j/i+i] X Z, Ui = Ib + i ■ (tb — Ib)/N. Then we verified that 

N 

IJ P{B,) C B. (30) 

We used a rigorous ODE solver of order 25 from the CAPD library which imple¬ 
ments the Hermite-Obreshkov algorithm proposed in [29]. Rigorous enclosure 
for P{B) returned by our routine is shown in Fig. 1101 

0.034 

0.033 

0.032 

z 0.031 

0.030 

0.029 

0.028 


Figure 10: The set B (in red) and a rigorous enclosure for P{B) (in yellow) 
obtained as the union of enclosures Ui=i — see (l30l) . 
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Chaos. Semiconjugacy of P^\-h to the Bernoulli shift is proved by means of 
the method of covering relations — the same as in [4^ but applied to different 
sets. It is sufficient to check the following geometric conditions 


T^yP‘^{y,z) 

< 

Im 

for(y, z) 

e 

{Im} X Z, 


> 

tn 

for(j/, z) 

e 

{tm} X Z, 

^yP^iy^z) 

< 

Im 

for(y, z) 

e 

{tn} X Z 

T^yP^iy^z) 

> 

Tn 

for(y, z) 

G 

{In} X Z, 


(31) 


where denotes the canonical projection onto the y coordinate. The geometry 
of these conditions is shown in Fig. [HI For the precise statement of a general 
theorem concerning, the method of covering we refer to |46] . 

The conditions m have been verified in direct computation. We did not 
need to subdivide any of the four edges of N and M that appear in (EH). 
Rigorous bounds on P'^{{Im} x Z), P'^{{rM} x Z), P'^{{lf^} x Z) and P^Htn} x 
Z), returned by our routine, are shown in Fig. 1111 



Figure 11: The sets M and N and rigorous enclosures of the images of their 
exit edges — see EH). 

Hyperbolicity and full conjugacy. Uniform hyperbolicity of Ti, is proved 
by means of the cone condition introduced in m- Here we use our algorithm 
for integration of variational equations. Derivatives with respect to initial con¬ 
ditions are necessary for computation of the derivative of Poincare map P^. Let 
Q be a diagonal matrix Q = Diag(A, y) with arbitrary coefficients satisfying 
A > 0 and y, < 0. It has been shown [IT] that if for all {y,z) G N Li M the 
matrix 

DP\y,zf ■Q-DP\y,z)-Q (32) 

is positive definite, then the maximal invariant set for P^ in NDM is uniformly 
hyperbolic. In our computations we used A = 1 and /i = —1000. 
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We uniformly subdivided both sets N and M onto 48 and 32 equal pieces, 
respectively (only y coordinate was subdivided). Then, each rectangle was sub¬ 
mitted to our routine that integrates the first order variational equations and 
computes derivative of the Poincare map P^. Given a rigorous bound of the 
derivative, we checked successfully the condition (15^ . Note that in the case of 
2x2 matrix it is easy to check positive definiteness of a matrix by the Sylvester 
criterion. I 


5.1 Comparison of time of computation. 

In the section, we discuss how the CPU-time needed for verification of the 
uniform hyperbolicity in Theorem |1] depends on the choice of the algorithm used 
to integrate variational equations. To this end, we did the following numerical 
experiment. For fixed parameters 

• TO — the order of numerical method, 

• tol — truncation error per one step of the numerical method, 

• Alg — the algorithm used to integrate variational equations (C^-Lohner 
or C^-HO algorithm) 

we compute the following three numbers 

• g]\[{m,tol), gM(jn,tol) — minimal natural numbers, such that using algo¬ 
rithm Alg, the method of order to with the tolerance tol we were able to 
check the cone condition (ED) subdividing uniformly the sets N, M onto 
gN and gM parts, respectively, 

• t{Alg) — CPU time of checking the cone condition on both sets N and 
M with the algorithm and parameters as above. 

Let us emphasize that the vector held of the Rossler system m contains only 
one nonlinear term. Hence, we have c/ = 1, and this is almost the worst linear 
case for the C^-HO method when the complexity is dominated by the matrix 
operations and the expected time savings from the C^-HO method are smaller 
— see analysis in Section [3] 

In Table [1] we present results from this experiment. We see that in each 
case the C^-HO algorithm is faster than the C^-Lohner algorithm. Higher com¬ 
putational complexity of the C^-HO algorithm is compensated by signihcantly 
smaller truncation error. Therefore, a routine that predicts the time step (the 
same routine was used in both cases) returns larger time steps for the C^-HO 
algorithm, and in consequence, the total computing time is smaller. We also 
notice that in some cases decreasing the tolerance increases the number of subdi¬ 
visions gN and gM needed to check the cone condition — see for instance the row 
with TO = 14 and tol = This is a consequence of many heuristics made 

in the implementation (for instance reorganization of doubleton representation 
after reaching some threshold values). These heuristics make the algorithm 
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Ci-Lohner 

Ci-HO 

t(LO) 

t(HO) 

m 

tol 

9m 

9n 

t(LO) 

9m 

9n 

t(HO) 

10 

10-10 

39 

33 

0.90 

48 

32 

0.82 

1.10 

10 

10-12 

38 

31 

1.29 

37 

31 

1.02 

1.26 

10 

10-14 

25 

31 

1.69 

23 

30 

1.20 

1.41 

10 

10-16 

25 

28 

2.25 

24 

25 

1.67 

1.35 

14 

10-10 

45 

52 

1.13 

48 

48 

0.84 

1.34 

14 

10-12 

42 

39 

1.22 

41 

36 

0.90 

1.35 

14 

10-14 

47 

33 

1.65 

49 

33 

1.21 

1.36 

14 

10-16 

36 

32 

1.70 

36 

31 

1.38 

1.23 

18 

10-10 

63 

77 

1.67 

62 

56 

mjm 

1.40 

18 

10-12 

48 

56 

1.41 

63 

49 

mm 

1.12 

18 

10-14 

44 

43 

1.76 

47 

38 

mm 

1.50 

18 

10-16 

40 

37 

1.91 

54 

34 

mm 

1.33 

22 

10-10 

151 

95 

3.36 

101 

67 

1.93 

1.74 

22 

10-12 

78 

78 

2.44 

59 

58 

1.81 

1.35 

22 

10-14 

52 

61 

2.05 

61 

49 

1.56 

1.32 

22 

10-16 

45 

41 

1.80 

47 

47 

1.53 

1.17 


Table 1: Comparison of C^-Lohner and C^-HO algorithms. 
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discontinuous with respect to parameters. Moreover, decreasing the tolerance 
increases number of time steps needed to compute a full trajectory. This may 
result in weaker control of unavoidable wrapping effect. 


6 Conclusions. 

Since the C^-Lohner algorithm appeared [47] it has been proved to be very useful 
in rigorous analysis of ODEs. In this paper, we proposed an efficient alternative 
for this algorithm and we provided free implementation of both C^-Lohner and 
C^-HO algorithms available as a module of the CAPD library |5]. Numerical 
tests show that the C^-HO algorithm is slightly faster than the widely used 
C^-Lohner algorithm. We have shown that the C^-HO algorithm may be faster 
in practical applications. This is not very important when the total time of 
computation is counted in seconds, as we have seen in Section |S| Any progress 
matters, however, if a problem requires hundreds or thousands CPU hours: for 
example verification of the existence of an uniformly hyperbolic attractor of the 
Smale-Williams type [41] or the coexistence of chaos and hyperchaos in the 4D 
Rossler system 13 Eg. In the computation reported in [3 Eg, the proposed 
C^-HO algorithm has been used. 

In [45] , an algorithm for integration of higher order variational equations is 
presented. Ideas from Section |3 can be directly used to design and implement 
an algorithm, let us call it C’'-HO, with the C’’-Lohner method as a predictor 
step. This requires encoding rather than theoretical effort and, we hope, this 
implementation will be available soon as part of the CAPD library [B]. 
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